#ifndef prop_H
#define prop_H

//#include <iostream>

#include "geometry.h"
#include "matrixes.h"
#include "mesh.h"
#include "boundary.h"
#include "exports.h"

namespace ABAN
{
    class EulerBackwardPropellant
    {
        public://protected but boundarycondition is friend
        FloatType nu,dt,density;

        SparseMatrix *velsCoeffMatrix[noOfAxis],*pCoeffMatrix;
        FloatVector *velsB[noOfAxis],*pB;

        DoubleDatumVector *dummyCellData;

        void createEqMatrixes(SeqType nn);
        void destroyEqMatrixes();

        virtual void applyBoundariesOnUMatrix(CellContainer* cells);
        virtual void applyBoundariesOnPSourceTerm(CellContainer* cells);
        virtual void applyBoundariesOnCorrectedU(CellContainer* cells, bool updateFromUHat);
        virtual void applyBoundariesOnPressureGrad(CellContainer* cells);

        virtual void calculatePressureGradients(CellContainer* cells);
        virtual void eqsGenerator(CellContainer* cells);
        virtual void generatePressureSourceTerm(CellContainer* cells);
        virtual void updateUS(CellContainer* cells);
        virtual void updateUN(CellContainer* cells);

        public:
        EulerBackwardPropellant(FloatType anu, FloatType aDensity, FloatType adt)
            :nu(anu), dt(adt), density(aDensity) {}
        virtual ~EulerBackwardPropellant(){}
        void calc(CellContainer *cells);
    };
    //----------------------------------------------------------------------------

    void EulerBackwardPropellant::calculatePressureGradients(CellContainer* cells)
    {
        for (CellIterator cellIter=cells->begin(); cellIter!=cells->end(); cellIter++)
        {
            Vector totalCorrection;
            SeqType pCellNo=cellIter->no;
            Vector rP=cellIter->center;
            for (SeqType face=dirXP; face<facesPerCell; face++)
            {
                FloatType coeff=(face<noOfAxis?1.0:-1.0);// cos(n<e)
                SeqType dir = (face<noOfAxis?face:(face-noOfAxis));
                Vector gP_P=cellIter->LSEgradP;
                for (CellNeighborhoodSetIterator neighbor=cellIter->neighbors[face]->begin();
                    neighbor!=cellIter->neighbors[face]->end();
                    neighbor++)
                {
                    SeqType nCellNo = (*neighbor)->no;
                    Vector rN = (*neighbor)->center;
                    Vector rF = cellIter->sharedFaceCenter((*neighbor),dir);
                    FloatType areaF = cellIter->sharedArea((*neighbor),dir);
                    FloatType addingVal=(rF[dir]-rP[dir])*(*cells->p)[nCellNo]+(rN[dir]-rF[dir])*(*cells->p)[pCellNo];
                    Vector gP_N = (*neighbor)->LSEgradP;

                    for(SeqType otherDir=dirXP;(otherDir<noOfAxis);otherDir++)
                        if (otherDir!=dir)
                            addingVal+=(rF[dir]-rP[dir])*(rF[otherDir]-rN[otherDir])*gP_N[otherDir] +
                                       (rN[dir]-rF[dir])*(rF[otherDir]-rP[otherDir])*gP_P[otherDir] ;

                    addingVal*=coeff*areaF/(rN[dir]-rP[dir]);
                    totalCorrection[dir]+=addingVal;
                }//neighb for
            }// face for
            (*cells->pGradients[dirXP])[pCellNo]=totalCorrection[dirXP];
            (*cells->pGradients[dirYP])[pCellNo]=totalCorrection[dirYP];
            (*cells->pGradients[dirZP])[pCellNo]=totalCorrection[dirZP];
        }// cell for
    }

    //----------------------------------------------------------------------------
    void EulerBackwardPropellant::updateUS(CellContainer* cells)
    {
        for (CellIterator cellIter=cells->begin(); cellIter!=cells->end(); cellIter++)
        {
            SeqType pCellNo=cellIter->no;
            (*cells->vels[dirXP])[pCellNo]+=(*cells->pGradients[dirXP])[pCellNo]*dt/density/cellIter->volume;
            (*cells->vels[dirYP])[pCellNo]+=(*cells->pGradients[dirYP])[pCellNo]*dt/density/cellIter->volume;
            (*cells->vels[dirZP])[pCellNo]+=(*cells->pGradients[dirZP])[pCellNo]*dt/density/cellIter->volume;
        }// cell for
    }
    //----------------------------------------------------------------------------
    void EulerBackwardPropellant::updateUN(CellContainer* cells)
    {
        for (CellIterator cellIter=cells->begin(); cellIter!=cells->end(); cellIter++)
        {
            SeqType pCellNo=cellIter->no;
            (*cells->velsN[dirXP])[pCellNo]=(*cells->vels[dirXP])[pCellNo]-(*cells->pGradients[dirXP])[pCellNo]*dt/density/cellIter->volume;
            (*cells->velsN[dirYP])[pCellNo]=(*cells->vels[dirYP])[pCellNo]-(*cells->pGradients[dirYP])[pCellNo]*dt/density/cellIter->volume;
            (*cells->velsN[dirZP])[pCellNo]=(*cells->vels[dirZP])[pCellNo]-(*cells->pGradients[dirZP])[pCellNo]*dt/density/cellIter->volume;
        }// cell for
    }
    //----------------------------------------------------------------------------
    void EulerBackwardPropellant::createEqMatrixes(SeqType nn)
    {
        pCoeffMatrix = new SparseMatrix(nn);
        pB = new FloatVector(nn);

        for (SeqType i=dirXP; i<noOfAxis; i++)
        {
            velsCoeffMatrix[i]=new SparseMatrix(nn);
            velsB[i]=new FloatVector(nn);
        }
    }
    //----------------------------------------------------------------------------
    void EulerBackwardPropellant::destroyEqMatrixes()
    {
        delete pCoeffMatrix;
        delete pB;

        for (SeqType i=dirXP; i<noOfAxis; i++)
        {
            delete velsCoeffMatrix[i];
            delete velsB[i];
        }
    }
    //----------------------------------------------------------------------------
    void EulerBackwardPropellant::calc(CellContainer *cells)
    {
        calculatePressureGradients(cells);
        applyBoundariesOnPressureGrad(cells);
        dummyCellData=new DoubleDatumVector(cells->size());

        for (SeqType i=0; i<300; i++)
        {
            char outName[80];
            createEqMatrixes(cells->size());

            cells->updateGradients();
            eqsGenerator(cells);
            applyBoundariesOnUMatrix(cells);

            //solving matrixes
            SparseMatrix::SORSolver(velsCoeffMatrix[0], velsB[0], cells->vels[dirXP], cells->velsN[dirXP],1.0,1e-7,100);
            SparseMatrix::SORSolver(velsCoeffMatrix[1], velsB[1], cells->vels[dirYP], cells->velsN[dirYP],1.0,1e-7,100);
            SparseMatrix::SORSolver(velsCoeffMatrix[2], velsB[2], cells->vels[dirZP], cells->velsN[dirZP],1.0,1e-7,100);

            updateUS(cells);
//            applyBoundariesOnCorrectedU(cells,true);

            cells->updateGradients();
            generatePressureSourceTerm(cells);
            applyBoundariesOnPSourceTerm(cells);


/*  printf("ooo  %f\n",(*pB)[0]);
  getchar();*/
            SparseMatrix::SORSolver(pCoeffMatrix, pB, cells->p, cells->p,1.7,1e-10,300);
//             SparseMatrix::BiCGStabSolver(pCoeffMatrix, pB, cells->p, cells->p);

//             getchar();

            if (!cells->hasDirichletPBC)
            {
                SeqType size=cells->p->size();
                for (SeqType row=0; row<size; row++)
                    (*cells->p)[row]-=(*cells->p)[size-1];
            }

            cells->updateGradients();
            calculatePressureGradients(cells);
            applyBoundariesOnPressureGrad(cells);

            updateUN(cells);
            for (SeqType row=0;row<dummyCellData->size();row++) (*dummyCellData)[row]=(*cells->velsN[0])[row];//must be deleted
//            applyBoundariesOnCorrectedU(cells,false);

            sprintf(outName,"./out/cav%d.vtk",(SeqType)((FloatType)i/10.0));
            if (i%10==0) VTKExport::run(cells,outName/*,dummyCellData*/);
            std::cout<<i<<std::endl<<std::flush;

            destroyEqMatrixes();
        }

        delete dummyCellData;
    }
    //----------------------------------------------------------------------------
    void EulerBackwardPropellant::generatePressureSourceTerm(CellContainer* cells)
    {
        for (CellIterator cellIter=cells->begin(); cellIter!=cells->end(); cellIter++)
        {
            FloatType totalFlux=0.0;
            SeqType pCellNo=cellIter->no;
            Vector rP=cellIter->center;
            for (SeqType face=dirXP; face<facesPerCell; face++)
            {
                FloatType coeff=(face<noOfAxis?1.0:-1.0);// cos(n<e)
                SeqType dir = (face<noOfAxis?face:(face-noOfAxis));
                Vector gU_P=cellIter->LSEgradV[dir];
                FloatType addingVal=0.0;

                SeqType nCellNo=-1;
                for (CellNeighborhoodSetIterator neighbor=cellIter->neighbors[face]->begin();
                    neighbor!=cellIter->neighbors[face]->end();
                    neighbor++)
                {
                    nCellNo = (*neighbor)->no;
                    Vector rN = (*neighbor)->center;
                    Vector rF = cellIter->sharedFaceCenter((*neighbor),dir);
                    FloatType areaF = cellIter->sharedArea((*neighbor),dir);
                    addingVal=(rF[dir]-rP[dir])*(*cells->vels[dir])[nCellNo]+(rN[dir]-rF[dir])*(*cells->vels[dir])[pCellNo];
                    Vector gU_N=(*neighbor)->LSEgradV[dir];

                    for(SeqType otherDir=dirXP;(otherDir<noOfAxis);otherDir++)
                        if (otherDir!=dir)
                        {
                            addingVal+=(rF[dir]-rP[dir])*(rF[otherDir]-rN[otherDir])*gU_N[otherDir] +
                                       (rN[dir]-rF[dir])*(rF[otherDir]-rP[otherDir])*gU_P[otherDir] ;

                        }
                    addingVal*=coeff*areaF/(rN[dir]-rP[dir]);
                    totalFlux+=addingVal;
                }//neighb for

//                 if (pCellNo==0)
//                 {
//                   printf("face : %d    Neighb:%d    flux : %f\n",face,nCellNo,addingVal);
//                   getchar();
//                 }

            }// face for
            (*pB)[pCellNo]+=totalFlux/dt;
        }// cell for
    }
    //----------------------------------------------------------------------------
    void EulerBackwardPropellant::eqsGenerator(CellContainer* cells)
    {
        for (CellIterator cellIter=cells->begin(); cellIter!=cells->end(); cellIter++)
        {
            SeqType pCellNo=cellIter->no;
            FloatType cellVol=cellIter->volume;
            Vector rP=cellIter->center;

            cellIter->checked=cells->currCheckedVal;
            Vector gradientsU_P[/*eqNum*/noOfAxis],gradientsU_N[/*eqNum*/noOfAxis];
            Vector gradientsP_P,gradientsP_N;

            gradientsU_P[/*eqNum*/0]=cellIter->LSEgradV[dirXP];
            gradientsU_P[/*eqNum*/1]=cellIter->LSEgradV[dirYP];
            gradientsU_P[/*eqNum*/2]=cellIter->LSEgradV[dirZP];
            gradientsP_P=cellIter->LSEgradP;

            velsCoeffMatrix[/*eqNum*/0]->setVal(pCellNo,pCellNo,cellVol/dt);
            velsCoeffMatrix[/*eqNum*/1]->setVal(pCellNo,pCellNo,cellVol/dt);
            velsCoeffMatrix[/*eqNum*/2]->setVal(pCellNo,pCellNo,cellVol/dt);
            pCoeffMatrix->setVal(pCellNo,pCellNo,0.0);

            (*velsB[/*eqNum*/0])[pCellNo]=(cellVol/dt)*(*cells->velsN[dirXP])[pCellNo]-(*cells->pGradients[dirXP])[pCellNo]/density;
            (*velsB[/*eqNum*/1])[pCellNo]=(cellVol/dt)*(*cells->velsN[dirYP])[pCellNo]-(*cells->pGradients[dirYP])[pCellNo]/density;
            (*velsB[/*eqNum*/2])[pCellNo]=(cellVol/dt)*(*cells->velsN[dirZP])[pCellNo]-(*cells->pGradients[dirZP])[pCellNo]/density;
            (*pB)[pCellNo]=0.0;

            for (SeqType face=dirXP; face<facesPerCell; face++)
                for (CellNeighborhoodSetIterator neighbor=cellIter->neighbors[face]->begin();
                    neighbor!=cellIter->neighbors[face]->end();
                    neighbor++)
                {
                    SeqType nCellNo = (*neighbor)->no;
                    Vector rN = (*neighbor)->center;
                    SeqType dir = (face<noOfAxis?face:(face-noOfAxis));
                    FloatType deltaRDir = (rN[dir]-rP[dir]);
                    Vector rF = cellIter->sharedFaceCenter((*neighbor),dir);
                    FloatType areaF = cellIter->sharedArea((*neighbor),dir);
                    FloatType coeff=(face<noOfAxis?1.0:-1.0);// cos(n<e)

                    gradientsU_N[/*eqNum*/0]=(*neighbor)->LSEgradV[dirXP];
                    gradientsU_N[/*eqNum*/1]=(*neighbor)->LSEgradV[dirYP];
                    gradientsU_N[/*eqNum*/2]=(*neighbor)->LSEgradV[dirZP];
                    gradientsP_N=(*neighbor)->LSEgradP;

                    //Normal to face Velocity (it is scalar - means that flow exits from face and + ...)

                    FloatType Uf=(rF[dir]-rP[dir])*(*cells->vels[dir])[nCellNo]+(rN[dir]-rF[dir])*(*cells->vels[dir])[pCellNo];
                    FloatType addingPart=((*cells->p)[nCellNo]-(*cells->p)[pCellNo]);
                    for(SeqType otherDir=dirXP;(otherDir<noOfAxis);otherDir++)
                        if (otherDir!=dir)
                        {
                            Uf+=(rF[dir]-rP[dir])*(rF[otherDir]-rN[otherDir])*gradientsU_N[dir][otherDir]+
                                (rN[dir]-rF[dir])*(rF[otherDir]-rP[otherDir])*gradientsU_P[dir][otherDir];
                            addingPart+=gradientsP_N[otherDir]*(rF[otherDir]-rN[otherDir])-
                                        gradientsP_P[otherDir]*(rF[otherDir]-rP[otherDir]);

                        }
                    Uf*=coeff/deltaRDir;
                    Uf-=addingPart/fabs(deltaRDir)/density*dt;
//Uf=0.0;//!!!



                    //Diffusive  active part diffActivePartForPCell,diffActivePartForNCell
                    FloatType diffActivePartForPCell,diffActivePartForNCell;
                    diffActivePartForNCell=-(diffActivePartForPCell=-nu*(1.0/fabs(deltaRDir))*areaF);
//diffActivePartForNCell=diffActivePartForPCell=0.0;
                    //Diffusive  active part for P
                    FloatType pressureDiffActivePartForPCell,pressureDiffActivePartForNCell;
                    pressureDiffActivePartForNCell=-(pressureDiffActivePartForPCell=-(1.0/density)*(1.0/fabs(deltaRDir))*areaF);
//pressureDiffActivePartForNCell=pressureDiffActivePartForPCell=0.0;
                    //Convective active part convActivePartForPCell,convActivePartForNCell
                    FloatType convActivePartForPCell,convActivePartForNCell;
                    convActivePartForPCell=(rN[dir]-rF[dir])/deltaRDir * Uf * areaF;
                    convActivePartForNCell=(rF[dir]-rP[dir])/deltaRDir * Uf * areaF;
//convActivePartForNCell=convActivePartForPCell=0.0;
                    //Putting active values in coefficient matrtixes
                    for (SeqType eqNum=dirXP;eqNum<noOfAxis;eqNum++)
                    {
                        velsCoeffMatrix[eqNum]->setVal(pCellNo,pCellNo,
                                                       velsCoeffMatrix[eqNum]->getVal(pCellNo,pCellNo) +
                                                       convActivePartForPCell-diffActivePartForPCell);
                        velsCoeffMatrix[eqNum]->setVal(pCellNo,nCellNo,
                                                       convActivePartForNCell-diffActivePartForNCell);
                    }
                    pCoeffMatrix->setVal(pCellNo,pCellNo,
                                         pCoeffMatrix->getVal(pCellNo,pCellNo)+pressureDiffActivePartForPCell);
                    pCoeffMatrix->setVal(pCellNo,nCellNo,
                                         pCoeffMatrix->getVal(pCellNo,nCellNo)+pressureDiffActivePartForNCell);

                    FloatType diffLaggedPartForPCell[noOfAxis]={0.0,0.0,0.0},
                              diffLaggedPartForNCell[noOfAxis]={0.0,0.0,0.0};
                    FloatType convLaggedPartForPCell[noOfAxis]={0.0,0.0,0.0},
                              convLaggedPartForNCell[noOfAxis]={0.0,0.0,0.0};
                    FloatType pressureDiffLaggedPartForPCell=0.0,
                              pressureDiffLaggedPartForNCell=0.0;

                    for (SeqType eqNum=dirXP;eqNum<noOfAxis;eqNum++)
                    {
                        for(SeqType otherDir=dirXP;(otherDir<noOfAxis);otherDir++)
                            if (otherDir!=dir)
                            {
                                //diffusive  lagged part diffLaggedPartForPCell,diffLaggedPartForNCell
                                diffLaggedPartForPCell[eqNum]+= nu*(rF[otherDir]-rP[otherDir])*areaF*
                                                                gradientsU_P[eqNum][otherDir] / deltaRDir;
                                diffLaggedPartForNCell[eqNum]+=-nu*(rF[otherDir]-rN[otherDir])*areaF*
                                                                gradientsU_N[eqNum][otherDir] /deltaRDir;


                                //Convective lagged part convLaggedPartForPCell,convLaggedPartForNCell
                                convLaggedPartForPCell[eqNum]+=(rN[dir]-rF[dir])*
                                                               (rF[otherDir]-rP[otherDir])*
                                                                areaF * Uf *
                                                                gradientsU_P[eqNum][otherDir] /
                                                                deltaRDir;
                                convLaggedPartForNCell[eqNum]+=(rF[dir]-rP[dir])*
                                                               (rF[otherDir]-rN[otherDir])*
                                                                areaF * Uf *
                                                                gradientsU_N[eqNum][otherDir] /
                                                                deltaRDir;
                            }

                        //Putting lagged values in b matrtixes
                        (*velsB[eqNum])[pCellNo]+=
                                                (diffLaggedPartForPCell[eqNum]+convLaggedPartForPCell[eqNum])+
                                                (diffLaggedPartForNCell[eqNum]+convLaggedPartForNCell[eqNum]);
                    }//of eqNum for

                    for(SeqType otherDir=dirXP;(otherDir<noOfAxis);otherDir++)
                        if (otherDir!=dir)
                        {
                            //diffusive  lagged part diffLaggedPartForPCell,diffLaggedPartForNCell
                            pressureDiffLaggedPartForPCell+=+1.0/density*(rF[otherDir]-rP[otherDir])*areaF*
                                                            gradientsP_P[otherDir] / fabs(deltaRDir);
                            pressureDiffLaggedPartForNCell+=-1.0/density*(rF[otherDir]-rN[otherDir])*areaF*
                                                            gradientsP_N[otherDir] / fabs(deltaRDir);
                        }
                    (*pB)[pCellNo]+=pressureDiffLaggedPartForPCell+pressureDiffLaggedPartForNCell;

                }//of face & neighb for
        }//of cellIter for

    }//of generate equations
    //----------------------------------------------------------------------------
    void EulerBackwardPropellant::applyBoundariesOnUMatrix(CellContainer* cells)
    {
        for (BoundaryContainer::iterator i=cells->boundary->begin();
             i!=cells->boundary->end();
             i++)
             (*i)->applyEffectOnUMatrixes(this,cells);
    }
    //----------------------------------------------------------------------------
    void EulerBackwardPropellant::applyBoundariesOnPSourceTerm(CellContainer* cells)
    {
        for (BoundaryContainer::iterator i=cells->boundary->begin();
             i!=cells->boundary->end();
             i++)
             (*i)->applyEffectOnPMatrixes(this,cells);
    }
    //----------------------------------------------------------------------------
    void EulerBackwardPropellant::applyBoundariesOnCorrectedU(CellContainer* cells, bool updateFromUHat)
    {
        for (BoundaryContainer::iterator i=cells->boundary->begin();
             i!=cells->boundary->end();
             i++)
             (*i)->applyEffectOnCorrectedU(this,cells,updateFromUHat);
    }
    //----------------------------------------------------------------------------
    void EulerBackwardPropellant::applyBoundariesOnPressureGrad(CellContainer* cells)
    {
        for (BoundaryContainer::iterator i=cells->boundary->begin();
             i!=cells->boundary->end();
             i++)
             (*i)->applyEffectOnPressureGrad(this,cells);
    }
    //----------------------------------------------------------------------------
}
#endif
